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Abstract. - A calculation of the rotational So(0) frequencies in high pressure solid para- 
hydrogen is performed. Convergence of the perturbative series at high density is demonstrated 
by the calculation of second and third order terms. The results of the theory are compared with 
the available experimental data to derive the density behaviour of structural parameters. In 
particular, a strong increase of the value of the lattice constant ratio c/a and of the internuclear 
distance is determined. Also a decrease of the anisotropic intermolecular potential is observed 
which is attributed to charge transfer effects. The structural parameters determined at the 
phase transition may be used to calculate quantum properties of the rotationally ordered phase. 



The phase diagrams of hydrogen and deuterium demonstrate amazing phenomena caused 
by the quantum nature of both translations and rotations of the molecules in these solids [1,2]. 
At low pressure, both systems have an hexagonal close packed (hep) structure, with almost 
freely rotating molecules [3] . A transition to a broken symmetry phase where the molecules are 
oriented, and the equivalence of the sites is lost (BSP or phase II) occurs for low temperature 
at 28 GPa for D2 and 110 GPa for H2. [4-6] An appropriate description of the BSP is still 
lacking. The correct quantum treatment of nuclei is indeed crucial for the description of 
the BSP transition, given the extremely large isotopic effect on the transition pressure. For 
this reason a description on the basis of a molecular interaction potential offers numerous 
advantages with respect to, for example, ab initio molecular dynamics [7] (which generally 
treats the nuclear motion as classical), quantum Monte Carlo techniques (limited to zero 
temperature), [8] or other more sophisticated computation techniques, still prohibitive for 
heavy computer time needed [9]. 

An effective potential model for the anisotropic molecular interaction in the solid has never 
been derived, especially because a reliable theory to derive solid spectroscopic data was lacking. 
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In this paper we present a theory able to calculate the roton frequencies in solid hep para- 
hydrogen up to very high pressure. From the comparison with experimental data, we obtain 
information on the anisotropic components of the hydrogen-hydrogen intermolccular potential, 
and we derive, as a function of pressure, fundamental structural parameters as the lattice 
constant ratio c/a and the internuclear distance. Here we give a further demonstration that 
effective two-body potentials are well suited to explain experimental results in high pressure 
systems and, as such, they are possibly able to explain the transition to the rotationally 
ordered phase. 

The triplet structure of the So(0) Raman transition for para- hydrogen is known since long 
time [10], and is considered a clear evidence of the rotationally disordered hep lattice. Its 
pressure evolution, including its interaction and hybridisation with the phonon, occurring at 
intermediate densities, has been measured up to the transition to the BSP [4-6], and has been 
the subject of a renewed interest also recently [11]. 

At low density a first-order perturbative theory is usually employed for the theoretical 
description of the triplet. In such a picture the first excited rotational state (J = 2) broadens 
in a band of collective states (rotons). Raman transitions from the ground state are allowed 
only to five states which correspond to q = rotons belonging to symmetry representations 
Ai g , E\ g and E^ g of the rotation group of the crystal. The E states are doubly degenerate and 
thus the Raman spectrum results in a triplet [3]. The comparison of the first order perturbation 
theory with experiment has not been satisfactory and the applicability of perturbative methods 
to the compressed solid has been questioned [11, 12]. In the following we demonstrate that, 
going beyond first order, the perturbation expansion is convergent and the numerical results 
are reliable up to the density of the BSP transition. 

We consider a rigid hep lattice of rigid molecules. The zeroth order eigenfunctions of the 
J = 2 Raman active (q = 0) roton modes are given by 



where i runs over the molecules and N is the number of unit cells, /j, can take the values 0, 
±1 and ±2 that correspond to the Ai g , Ei ff , E2 9 representations respectively. The spherical 
harmonics Yhm are calculated with respect to a reference system having the z-axis coincident 
with the c-axis of the hep crystal. The anisotropic intermolecular potential is assumed to be 
pair-additive, and is expanded in spherical components in the usual way (see eqs. 10 and 11 
of ref. [13]). We consider here those spherical components that have been evaluated in the 
literature, [14-16] namely V202 (= V022) that takes into account crystal field effects, and V22L 
(with L = 0, 2, 4) that is responsible for rotonic states. 
First order energies for the excited state are written [3] 



where C{L\ L2 L; mi, m,2, m) are the Clebsch-Gordan coefficients, and 5; 1 ; 2 l are the lattice 
sums 
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The first order correction removes the degeneracy among states with different |/x|, giving rise 
to three Raman active transitions of frequencies i/| M |,with uq > vi > v\. The splitting is 
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symmetric (yo — v 2 = v 2 — v\ ) only if calculated at first order and neglecting all but the V224 
potential components. 

The general expression for second order corrections to the J — 2 state is [17] 



4?> = E K TT' -E^m < 4 > 



where s indicate any rotational state other than the single J — 2 excitation and $(s) and i? s 

are the zeroth order wave- functions and energies for the rotational state s. The second equality 
( 21 

is a definition of ; (s). We have also calculated the expressions for the third order corrections 
both to the ground state and to the J = 2 rotational state. The calculation involves now two 
intermediate rotational states s\, s 2 [17] and the different terms are indicated with (si; S2) 
accordingly. 

Going beyond first order one finds that also the energy of the ground state changes, with a 
correction that is proportional to the number of molecules 2N, having therefore an extensive 
character. Anyway, both the second and third order corrections to the energy of the excited 
state contain quantities proportional to 2N, that are identical to the corresponding corrections 

of the ground state and thus cancel out when calculating the transition frequencies v^. For 

(2) 

the frequencies at second order the terms that contribute are 

hvf = c «'(22) + C W'(222) + e« (42) + e« (40) (5) 

where the prime means that the extensive contributions have been subtracted and the ar- 
gument of e here explicitly indicates the number and nature of rotational excitations of the 
intermediate state s (e.g. 222 indicates 3 rotational excitations with J = 2). The terms 
appearing in Eq. depend on products of the lattice sums 6*202 and S22L and on the new 
lattice sums Tll'm [18] 

Tlvm = 4tt^ V 22L (R tj )V 22L ,{R l3 )Y LM {n l3 )Yl, M {n l3 ). (6) 

At third order the terms that contribute are 

hv^ = e| 1 3) '(22;22) + e( l 3) '(222;222) + e(f ) '(42;42) + 

U 2) '(22) | e ' 2) '(222) | e < 2) (42) " 
" A 20 2A 20 A 40 

+ e( 1 3) '(222;22) + e| 1 3) '(222;42) + e( I 3) '(22;42) (7) 

where A40 and A20 are zero order energy differences of rotational states. The calculation 
requires the computation of double and triple lattice sums [18]. For the third order correction 
only (Eq. Q| we have neglected crystal field effects ^202- This is justified by the smallness of 
third order effects (see below). 

By the use of the analytical expressions for V224 and V202 given by [14] and the interpolation 
of the table given in [15] for V220 and V222 we have calculated the frequencies of the So(0) 
triplet, and analysed independently the first, second and third order correction terms. 

Even if the first order term exceeds 100 cm -1 at high density, the perturbative expansion 
is convergent up to 340 mol/1 and more, since the third order correction is always less than 
2 % of the unperturbed value v gas — 354.4 cm -1 . Writing explicitly the corrections in a 
typical high density case (p — 300 mol/1, c/a = ^/8/3), we have vq = v gas +94+14 + 5, 
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v \ = v gas — 36 + 43 — 2 and ^2 = +29 + 34+ 1 cm -1 . We have found that, even though 
the main contribution arises from the V224 potential component, which comprises the electric 
quadrupole-quadrupole (EQQ) interaction, the effect of the other potential terms are relevant 
(~ 20 cm -1 ) if c/a deviates a few percent from the ideal value. The effect of the crystal field 
term, V202 = V022 is negligible (only 0.02 %) for c/a — y/8/3. At high density, varying c/a of 
3 % may produce a change of 30 cm -1 for the absolute value of the frequencies and of about 
44 % of the splitting. The c/a ratio is therefore an important parameter to be considered 
when trying to reproduce the experimental rotational frequencies. 

It is important to notice an increase of the calculated average energy of the triplet at high 
density (this is true even considering only first order term if V220 is included in the calculation) 
and an asymmetric reduction of the splitting. Results of the calculation, including the third 
order, with fixed rotational constant Bq — 59.246 cm -1 are reported in fig. ^ We have 
compared our theoretical results with experimental data [4-6,11], also reported in fig. ^ as 
a function of density p calculated from the reported pressure by using the room temperature 
equation of state, neglecting thermal pressure effects [19]. The experimentally observed overall 
decrease of the rotational frequencies above a certain density contrasts with the increase of 
the theoretical values and addresses for a decrease of Bq with density. It is also evident that 
above 150 mol/1 the calculated splitting overestimates the experimental one. This problem 
was already evident to other authors who have used only first order expansion and only EQQ 
interaction [11, 12]. Such a discrepancy in one case has been ascribed to the believed failure 
of the perturbative expansion above 40 GPa, and it is corrected below 40 GPa assuming a 
screening of the EQQ interaction in the solid, proportional to the square root of the refractive 
index [12]. By contrast, other authors assume that the first order perturbation theory is 
sufficient up to high pressure, and introduce a new crystal field term, that they calculate 
at each pressure ascribing to this one the whole discrepancy [11]. Both approaches are not 
completely convincing, because they are limited by the first order perturbative theory. We 
are in a much better position, since, due to the convergence of the perturbative series, we 
can confidently push our analysis to the whole experimentally investigated density range, 
considering the absolute value of the individual roton frequencies and not only the splitting. 

Our conclusion is that, in addition to Bq, also the anisotropic interaction energy has to be 
assumed to change significantly with increasing density in the solid phase. We have modelled 
this change in an effective way, by multiplying all the potential components by the same 
factor, (3. We have pushed our analysis even further, taking full advantage of the information 
contained in the experimental data of the triplet frequencies, to derive, at any density, the 
third unknown parameter entering in the calculation, namely, the lattice constant ratio c/a. 
The results of the analysis are reported in fig. ^ Letting only B to vary (c/a = y/8/3, 
f3 = 1 ) we cannot find an agreement for the individual frequencies, since the splitting is still 
overestimated (dashed line), while we can perfectly reproduce the data by adjusting B , (3, 
and c/a. The analysis has been performed using the expansion up to second order (solid line 
with dots), and checking consistency summing up the third order contribution (open circles). 

The density behaviour of (3 and c/a are reported in fig. Even though an undoubtedly 
large effect ((3 ~ 0.3 at high density), the attenuation of the pair anisotropic potential is not 
surprising. Moreover, the 70 % variation we have found is, in absolute value, much smaller than 
the difference between the gas [15] and the solid [24] isotropic potential. Effective anisotropic 
interaction has been occasionally used to calculate properties of high pressure hydrogen, and 
in particular the pressure where the transition to the BSP is expected. Using the gas value of 
the EQQ interaction has always produced too low values for the transition pressure, indicating 
an overestimation of the anisotropic interaction [25,26]. For HD, where the transition line 
has a particular reentrant behaviour [27], a correct estimation of the transition pressure can 
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Fig. 1 - Theoretical results with fixed parameters (solid lines) and comparison of the results of the 
calculation with the experiment. Solid grey squares and dots are the experimental data for hydrogen 
of ref.s [4, 11]. The dashed line and solid dots joined by a line represent the values calculated with 
different free parameters (see text). The empty circles are the values calculated including third order 
terms, with parameters corresponding to the solid dots. 



be obtained by a mean-field theory and considering the effect of a crystal field term that may 
be physically produced by a lattice distortion [28] . In fig. we report also the attenuation 
factor assumed by Loubeyre et al. [12] and that derived by Runge et al. [20] and used in a 
Path- Integral Monte Carlo calculation [29] . Both of them are of the same order of magnitude 
as ours. 

For what concerns the lattice constant ratio c/a, reported in fig. |2t>, we notice that, up 
to p ~200 mol/1 the spectroscopically determined value differs only slightly from either -Jo/3 
or the room temperature x-rays value, while, above, a marked increase is evident. A recent 
self-consistent analysis of lattice distortion concludes, on the basis of a mean-field approach, 
that the hep lattice remains ideal up to the BSP transition pressure, where a large variation of 
c/a occurs [30,31]. This increase occurs for different densities in the cases of H2 and D2 [18], 
thus demonstrating that it represents a precursor of the phase transition, rather than a mere 
density effect. At room temperature, indeed, this increase is not present [19], while an increase 
of the lattice constant c passing from the hep to the BSP phase in D2 at constant pressure 
has been noticed in a x-ray diffraction measurement [32]. 
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Fig. 2 Fig. 3 

Fig. 2 - (a) Attenuation factor /3 versus density. The values used in ref. [12] (triangles) and in ref. [20] 
(circles) are reported, (b) Lattice constants ratio c/a. The grey symbols are the room temperature 
x-rays determination [19]. 

Fig. 3 - Bond length for H2 as a function of density. The grey symbols are values determined in the 
gas at room temperature (dots ref. [21], triangles ref. [22] and squares ref. [23]); diamonds are data 
derived in ref. [12]. The big solid dot is the value derived from data in ref. [10] and the horizontal 
grey lines represent values derived from the rotational constants in the low density gas. 



In fig. we report the effective internuclear distance r rn = (r~ 2 ) 12 = {^tt^cBq/K)^ 1 / 2 . 
In the same figure, we show data referring to the gas at room temperature [21-23], to the iso- 
lated molecule [33] and to the low pressure solid [10]. The marked decrease of the internuclear 
distance, already seen in ref. [12], proceeds only up to about 220 mol/1, where a rapid increase 
takes place. At the BSP transition, it becomes even larger than the isolated molecule value, 
probably favoured by the onset of the orientational ordering. Our results show that the rise of 
r m and c/a are two important features to be interpreted as precursors of the BSP transition. 
The increase of r m (with an expected increase of the molecular quadrupole moment [3] ) and 
the decrease of the anisotropic potential seem contradictory due to the eminently quadrupolar 
character of the interaction responsible for the Srj(0) structure. Such a contradiction calls for 
the consideration of charge transfer effects which are present even at low pressure and are 
especially important for determining the anisotropic intermolecular interaction. The struc- 
tural parameters and interaction potential at the transition should enable one to calculate 
equilibrium structure and dynamical effects in the BSP. 



t. (JRAZZI et al.: 1 HE S> (UJ TRIPLET IN PRESSURIZED HYDROGEN 



I 



REFERENCES 

[I] I. F. Silvera, Rev. Mod. Phys. 52, 393 (1980). 

[2] H. K. Mao, and R. J. Hemley, Rev. Mod. Phys. 66, 671 (1994). 
[3] J. Van Kranendonk, Solid Hydrogen, Plenum Press, (1983). 
[4] I. F. Silvera and R. J. Wijngaarden, Phys. Rev. Lett. 47, 39 (1981). 
[5] R. J. Wijngaarden, PHD thesis, Amsterdam University (1982). 

[6] H. E. Lorenzana, I. F. Silvera, and K. A. Goettel, Phys. Rev. Lett. 64 1939, (1990). 

[7] J. Kohanoff, Phys. Rev. Lett. 78, 2783, (1997). 

[8] D. M. Ceperley and B. J. Alder, Phys. Rev. B 36, 2092 (1987). 

[9] H. Kitamura, S. Tsuneyuki, T. Ogitsu, and T. Miyake, Nature, 404, 259 (2000). 

[10] S. Bhatnagar, E. J. Allin, and H. Welsh, Can. J. Phys. 40, 9 (1962). 

[II] A. F. Goncharov, M. K. Strzhemechny, H. K. Mao, and R. Hemley, Phys. Rev. B 63, 064304 
(2001). 

[12] P. Loubeyre, M. Jean-Louis, and I. F. Silvera, Phys. Rev. B43, 10191 (1991). 

[13] F. Grazzi, M. Santoro, M. Moraldi, and L. Ulivi, Phys. Rev. B 66, 144303 (2002). 

[14] M. J. Norman, R. O. Watts, U. Buck, J. Chem. Phys., 81, 3500 (1984), where the reported 
values of C5 and CV for V224 are exchanged, to reproduces correctly the long-range behavior of the 
potential, without altering the conclusions presented in the present reference (U. Buck, private 
communication) . 

[15] J. Schaefer, W. E. Koehler, Z. Phys. D, 13, 217 (1989). 

[16] P. Diep, J. K. Johnson, J. Chem. Phys., 112, 4465 (2000). 

[17] L. Landau E. Lifchitz, Mecantque Quanttque, Mir Editions, (1974). 

[18] F. Grazzi, M. Moraldi, and L. Ulivi, to be published. 

[19] P. Loubeyre, R. LeToullec, D. Hausermann, M. Hanfland, R. J. Hemley, H. K. Mao, and L. W. 

Finger, Nature, 383, 702 (1996). 
[20] K. J. Runge, M. P. Surh, C. Mailhiot, and E. L. Pollock Phys. Rev. Lett. 69, 3527 (1992). 
[21] A. D. May, V. Degen, J. C. Stryland, and H. L. Welsh, Can. J. Phys. 39, 1769 (1961). 
[22] N. E. Moulton, G. H. Watson, W. B. Daniels, and D. M. Brown, Phys. Rev.A37, 2475 (1988). 
[23] L. Ulivi, M. Zoppi, L. Gioe, and G. Pratesi, Phys. Rev. B58, 2383 (1998). 
[24] I. F. Silvera, V. V. Goldman, J. Chem. Phys. 69, 4209 (1978). 
[25] W. England, J. C. Raich, and R. D. Etters, J. Low Temp. Phys. 22, 213 (1976). 
[26] W. B. J. M. Janssen and A. van der Avoird Phys. Rev. B 42, 838 (1990). 
[27] F. Moshary, N. H. Chen, and I. F. Silvera, Phys. Rev. B48, 12613 (1993). 

[28] Yu. A. Freiman, S. M. Tretyak, A. Jezowski, and R. J. Hemley J. Low Temp. Phys. 113, 723 
(1998). 

[29] T. Cui, E. Cheng, B. J. Alder, and K. B. Whaley Phys. Rev. B 55, 12253 (1997). 

[30] Yu. A. Freiman, S. M. Tretyak, A. Jezowski, and R. J. Hemley J. Low Temp. Phys. 122, 537 

(2001) . 

[31] Yu. A. Freiman, S. M. Tretyak, A. Jezowski, and R. J. Hemley J. Low Temp. Phys. 126, 703 

(2002) . 

[32] P. Loubeyre, private communication (2001). 

[33] J. W. Foltz, D. H. Rank, and T. A. Wiggins, J. Mol. Spectrosc. 21, 203 (1966). 



